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Densities of Levy walks and the corresponding fractional equations 


Marcin Magdziarz Tomasz Zorawik ^ 


Abstract 

In this paper we derive explicit formulas for the densities of Levy walks. Our results cover both 
jump-first and wait-first scenarios. The obtained densities solve certain fractional differential 
equations involving fractional material derivative operators. In the particular case, when the 
stability index is rational, the densities can be represented as an integral of Meijer G function. 
This allows to efficiently evaluate them numerically. Our results show perfect agreement with 
the Monte Carlo simulations. 


1 Introduction 

The Continuous Time Random Walk (CTRW) is a stochastic process determined uniquely by i.i.d. 
random variables Ji, J 2 , • • • representing consecutive jumps of the random walker, and i.i.d. positive 
random variables Ti,T 2 ,... representing waiting times between jumps, m- The trajectories of 
CTRW are step functions with intervals Tj and jumps Jp 

Levy walk is a particular case of CTRW satisfying the following additional condition | Jj| = Tj for 
every i G N (length of the jump equal to the length of the preceding waiting time). Levy walks were 
dehned for the hrst time in ISHEI] in the framework of generalized master equations. Since then, 
they became one of the most popular models in statistical physics with large number of important 
applications. Levy walks have been found to be excellent models in the description of various real-life 
phenomena and complex anomalous systems. The most striking examples are: transport of light in 
optical materials [T], foraging patterns of animals mmm. epidemic spreading HE], human travel 
ElIH], blinking nano-crystals |16) . and fluid flow in a rotating annulus m- 

The main idea underlying the Levy walk is the spatial-temporal coupling, which is manifested 
by the condition | Jj| = Tj. Thus, even if we assume that the distribution of jumps Yi is heavy-tailed 
with diverging moments, the Levy walk itself has hnite moments of all orders (intuitively, long jumps 
are penalized by requiring more time to be performed). This is very different from a-stable Levy 
processes with a < 2, which have inhnite second moment. These desirable properties make Levy 
walk particularly attractive for physical applications. 

Let us now recall the formal dehnition of Levy walk. Let Tj, i = 1, 2,..., be the sequence of i.i.d. 
positive random variables, representing the waiting times of the walker. Assume that they belong 
to the domain of attraction of one-sided a-stable law, i.e. 

[nt] 

^ SUt) 

i=l 

as n ^ 00 . Here, Sait) is the a-stable subordinator with the Laplace transform j9] 

IE(exp{—sSQ,(t)}) = exp{—ts"}, 0 < a < 1. 
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Denote by 


k 

N{t) = max{A: > 0 Ti < t} 

i=l 

the corresponding counting process. Next, define the jumps of the walker as 


T — T T 

•Ji — 


where are i.i.d. random variables, which are assumed independent of the sequence of 

waiting times Tj. Each /j governs the direction of the jump, i.e. 

P(/i = 1) = p , P(/i = -1) = l- p, 


0 < p < 1. Thus, the walker jumps up with probability p and down with probability 1 — p. Clearly, 
the condition | Jj| = Tj is satisfied. The sequence of jumps Ji defined above belongs to the domain 
of attraction of a-stable distribution |13) 

[nt] 

i=l 

as n ^ oo. The process La{t) is the a-stable Levy motion with the corresponding Fourier transform 
E(exp{iA:LQ,(t)}) 

= exp{—cos(7ra/2)(l — i{2p — 1) tan{TTa/2)sgn{k))}. 


Later on we will also use the corresponding left-limit process 

limL„(s). 

s/'t 


Finally, the process 


N{t) 

m = 


2=1 

is called Levy walk. R{t) is also known as wait-first Levy walk in the literature |25) . since the walker 
at the beginning of its motion (t = 0) first waits and then performs the jump. 

As shown in |13) . R{t) obeys the following scaling limit in distribution 


AM ^ 


( 1 . 2 ) 


as n —)• oo. Here X{t) is the right-continuous version of the process L~{S~^{t)). Moreover, S~^(t) 
is the inverse of Sa{t), i.e. S^^{t) = inf{r > 0 : Sa{T) > t}. Additionally, the instants of jumps as 
well as the respective jump lengths of the processes La{t) and Sa{t) are exactly the same. Also, the 
probability density function (PDF) pt{x) of X{t) satisfies the following fractional equation jlOl 113) 


p 




+ (1 -P) 


\ / 


Pt{x) = (5o(x) 


r(l - a)' 


(1.3) 


Here, the operators =p are the fractional material derivatives introduced in |23) . In the 

Fourier-Laplace space they are given by 


XyCt 


\dt^^) 



{s^ikYf{k,s). 
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In what follows, we will find the explicit solution of (|1.3p by determining the PDF pt{x) of the limit 
process X{t) given in (II.2p . 

We will also consider the so-called jump-first Levy walk m 

N(t)-ei 

R{t)= 

i=l 

It has the following scaling limit in distribution ddS] 


R{nt) d 


Y{t) = L^{S-\t)) 


n 


(1.4) 


as n 


oo. Moreover, the PDF wt{y) of Y{t) satishes the fractional equation of the form [lOI 114) 


p 


^ _ A 

dt dy 


+ (1 -p) 


dt dy 


wtiy) = 


a 


r(l - a) 


6o{y — u)u “ 


(1.5) 


In the next two sections we will derive explicit formulas for the PDFs of the limit processes X{t) 
and Y(t) , respectively. This way we will obtain solutions of fractional equations (II.3p and (11.51) . 


2 Densities of wait-first Levy walks 

Let us start with the wait-hrst scenario. In the result below we determine the PDF of the limit 
processes X{t) from eq. (11.21) . 

Theorem 2.1 //p G (0; 1) a G (0,1) then the PDF pt{x) of the process X{t) = L~{S~^{t)) 
equals 


ptix) = —i 


1 


a 


2pa (1 — p)a r(l a) 

j-t yoo ^ 1 -a (W -\- x 

J\x\ Jo 


a 1 1 Z I Ta ( dzdwl(^o^t){\x\), 


( 2 . 1 ) 


'\x\Jo ^) 2pa J \2(1 — p) “ 

where ra{x) is a density of a positive a-stable random variable with the Laplace transform 

= e"““. 

Proof. We start with reminding the following equality of distributions (see |13)1: 
{La{t),Sa{t)) 

= (p^S^^\t) - (1 -p)-5®(t),p«5i^)(t) -h (1 -p)“^ 


( 2 . 2 ) 


a \^J 


where Sa\t) and Sa\t) are independent copies of Sa{t). Using this fact it was shown in |13) that 
the Levy measure i^[La,Sa) {La{t)-,Sa{t)) equals 


v{dx,ds) = p5s{dx)vsc{ds) + (1 - p)5-s{dx)vSoJds), 


(2.3) 
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where ^Sa{ds) = r(i-Q) ^ ^ ^ds. Hence, an infinitesimal generator A of the process {La{t)j Sa{t)) 
equals 


^f{x,s)= / [f{x + y,s + w) - f{x,s)]iy(^L^ Sc)idy,dw) 

7r2 

= / [f{x + y,s + w) - f{x,s)] 

{vdw{dy)vsc.{dw) + (1 - p)6-w{dy)vsc.{dw)) 


(2.4) 


InCZ] the authors consider (among other processes) the 2-dimensional process {X{t), V(t—)), where 
V (t) is the age process counting time that passed since the last jump of X{t). Equation (12.41) together 
with Remark 4.2 in the mentioned paper provides us with the joint distribution of {X{t),V{t—)) 

P {X{t) = dx,V{t-) = dv) = X [v,oo))U{dx,t - dv)l{o<^<t}: (2-5) 

where U{dx,ds) is the 0-potential measure of the process {La(t), Sa{t)), defined as 

poo 

U{dx,ds)= / P {La{u) = dx, Sa{u) = ds) du. 

Jo 

The process V{t) is beyond our interest, but we will later obtain the distribution of X{t) as a marginal 
distribution of {X{t),V{t—)). This is the reason why we turn now our attention to Ea. (|2.5p . We 
have 


X [u,oo)) =pus^{[v,oo)) + {l-p)us^{[v,oo)) = us^{[v,oo)) 




( 2 . 6 ) 


Jjj r(i — a) r(i — a) 

Furthermore, a density u{x, s) of the potential measure U can be calculated as 

poo 

u{x,s)= / Wu{x,s)du, 

Jo 

where wt{x,s) is the PDF of the process {La{t), Sa{t)). From Eq. (12.2p we get 
P (La(t) = (ix,5a(t) = ds) = P 


(2.7) 


2pc 


2(1-p)' 


The Jacobian determinant in the above formula for the linear transformation of {dx, ds) equals 
2 pi/a(i-p)i/a • Taking into account the independence of Sa\t) and Sa\t) we express wt{x,s) in 
terms of qt{s) - the PDF of Sait). 


1 


Wt{x,s) = -- 

2pa (1 — p) 


s + x\ 

rQt —^ qt 


s — X 


2pa J V2(l — p)< 
Moreover the self-similarity of Sa{t) provides us with the equation 


qu{x) = -^qi ( 


Ua \U°‘ 


1 


Ua 


Ua 


4 











here ra{x) is the density of the random variable 5 q( 1). We take into account the two above equations 
in Eg. (12.71) and substitute u = z~°^ to calculate the integral: 


m(x, s) = 


f 


L 

L 


0 2pa{\ — p^a 


1 ( S + X 

TT 1 ~ 

2pa [1 — p)a y 2p°‘ ^ 

1 ^ 

-u - ra 


Qu 


s — X 


2{l-p)^ 


du 


a 


0 2pa(\ — p)a 




s + X ^ \ 

- —U o‘ To 

( S + X 


s — X ^ \ , 
—u du 


.2(1-p)- 

s — X 


( 2 . 8 ) 


-z rn 


2pc 


.2(1 -p)“ 


—z dz. 


Finally we combine eqs. (1^ . (1^ . (1^ and integrate with respect to dv to obtain 

a 1 


Pt{x) = 


r(l a) 2pa (1 — p) 

pt—\x\ poo 

Jo Jo 




t — V + X \ ( t — V — X \ 

-j—z Tq - -z dzdv 

^ 2 p« J V 2 (l-P)“ J 


a 


r(l - «) 2p^{l-p)ii 

j-t rcxi ^1-a (W + X 

J\x\ Jo 


{t — wY 


Z Tn 


2pc 


w — X 


.2(1-P) 


-z dzdw 


We used here the fact that ra{x) vanishes outside the positive half-line and substituted v = t — zd 

One should mention here that the PDF of wait-first Levy walk in the extreme cases p = 0 and p = 1 
was already derived in m- 

For the special case a = ^ we get the following simple expression for pt{x): 

Corollary 2.1 When cx = ^ and p G (0,1) the PDF pt{x) can be expressed as 


\pl 

pt[x) = - p{l -p) 

TT 


(t- |x|)2 


-rl{0,i)(kl) 

{2p^t -I- (1 — 2p){t + x)) {2p^\x\ + (1 — 2p){x + |x|)) 2 

Proof. We have (see m) 

Taking into account this formula and using the property of Gamma function r(3/2) = ^r(l/2) = 
we can calculate the following integral: 


f 


1 

Z'2'i 


W + X 


w — X 


^1/2 

-3 

1 ( {w + x){w — x)\ 


z I dz 


dvr \ 4p2(l — pY 


f 


-5 

z 2 exp 


p‘^{w — x) -|- (1 — pY{w -|- x) 1 


2{w + x){w — x) 


(2.9) 


dz 


25 


p^{i — pY 


{p‘^{w — x) -|- (1 — pY{w + x)) 2 
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Substituting Eq. (|2.9p into Ea. (|2.1l) we arrive at 


Pt{x) = -P) 

V 27r 


f {t-w) 2 - pf + p^) w + {{I - pf - p"^) x) 2 du;l(o,i)(|x|). 

J\x\ 

One can notice that (below a, b and c are constants) 


d I —2 {c — w)2 
\(a + 6 c)(a + 6 tc )2 


= {c — w) 2 (a + hw) 2 , 


which finally implies 


pt{x) = - p{l-p) 

vr 




{2p^t + (1 — 2 p)(t + x)) {2p^\x\ + (1 — 2p){x + lx |))2 


l(o,t)(kl)- 


□ 


It is worth to mention here, that in a special case a = ^ and p = ^ we recover the known result 

m) 

Pt{x) = ^ [ it-w)-^w-^dwl^ot){\x\) = ^ 

J\x\ TTt\x\2 

Figures [T] and | 2 ] present the densities obtained from the above Corollary. 



X 


Figure 1: Plot of PDF pt{x) of the process (S'^ ^{t)) for different values of t, a = 0.5, p = 0.1. 

The drawback of using Theorem 12.11 to calculate pt{x) in the general case is the necessity of 
knowing values of ra{x). One can use some known algorithms to approximate ra{x), however results 
are not perfect. Computations are time-consuming and inaccurate. They differ from those obtained 
via Monte-Carlo methods. However, as the next Corollary shows, we can express pt(x) in the form 
of an integral from Meijer G function (see |19)b This special function is implemented in most of 
numerical packages, including Mathematica and Matlab. This representation is valid for rational a 
and can be used to approximate irrational a with a rational one. 
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Figure 2: For a = 0.5 we compare the obtained densities pi(x) for different values of p (red solid 
lines) with densities estimated using Monte Carlo methods (blue pluses). 


Corollary 2.2 If p G (0,1) and a = ^, where l,k then the PDF pt{x) of the process X(t) = 
equals 


Pt{x) 


2l a k+l Q, rt J 1 

pj^k-l [Cf Y{1 — a) {t — w)^ w — x 


p,k,k ( ( P I' w + x V -A(fc,fc(f-l)),A0,O) 
«+A:,/+A: ^ ytC _ X/ A(fc,0),-A(i,i(f-l)) 




where 


A(fc,a) 

A(/,6) 


is the Meijer G function (see m) and A{k, a) = 


a g+l 
A: ’ k ' 


a + 


a+2 


g+fe—1 


is a special list of k elements. 


Proof. For a = l/k where k, I are positive integers with k > I we can express r^ix) in terms of the 
Meijer G function (see |18)1: 


ri/k{x) 


^fikfi f}!_ -I 
(27r)(*^“*)/2 X \k^ 


A(i,o) \ 

A(fc,0) J- 
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This gives us 



( 27 r )(^“0 [w — x){w + x) Jq 


k 


4(1 — p)c‘pa 


(27r)(^“^) {w — x){w + x) Jq 



A{1,0) 

A{k,o) 


A(l,0) 

A(k,0) 


dz 


A(l,0) 

A(k,0) 


k 


A /-I \ \ 

A(1 — p)apa ( I \ W + X \ 


A{1,0) 

A{k,0) 


du 


(27r)(^“^) {w — x)(i(; + x) \k^ 

p 


G 


k^k I 

l+k,l+k I 


1 — p 


- / I \' 

“ w + x 


w — x 


2pa J 

-A{k,k{f-l)),A(lfi) 
A(fc,0)-A(i,i(f-l)) 


We substituted in the above equations z~^ = u. Moreover we used properties of Meijer G functions 
(see |19) ) to calculate the integral from multiplication of these functions. To end the proof we 
substitute the calculated above integral into Eq. m- n 


Figure |3] presents densities calculated here. Another method for calculating r^, where is not neces¬ 
sarily rational, was presented in |20) . 



Figure 3: For p = 0.25 we compare the densities pi{x) obtained in Corollary 12 .2 1 for different values 
of a (red solid lines) with densities estimated using Monte Carlo methods (blue pluses). 


3 Densities of first-jump Levy walks 

This section is devoted to the first-jump scenario. Below we find the PDF of the process Y{t) from 
m- At the same time this PDF solves the fractional equation (jl.Sp . 
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Theorem 3.1 If p € (0,1) and a G (0,1) then the PDF wt{y) of the process Y(t) = La{S^^{t)) 
equals: 


(i) if \y\ < t, then 

wt{y) = 


2^1-G 


a p 


r f 

J—t J\x\\/(t—v+x) J 0 


2(1 — p) « r(l — a) j—t J\x\\/{t—y+x) J 0 


W + X 


w — X 


-Z Tn 


—z\{y — x) ^ °‘dzdwdx 


0^(1 — p) 
+ —1 - - 


\ 2pa 
1-i rt rt 


ff 

J V J \x 


2paY{\ — Oi) Jy J\x\y{t+y—x) JQ 
( w + X \ 


2{l-p)^ 

poo 

/ 

Jo 


w — X 


-Z Tn 


—z\{x — y) ^ °‘dzdwdx, 


2pc 


2{l-p)^ 


(a) if y >t, then 

wt{y) = 


2 ^ 1 -- 


a p 


rt pt poo 


2(1 — p)c«r(l — a) J-tJ\x\Jo 

I w + X 


2pc 


-z rn 


^1 — OL 


w — X 


2{l-p)^ 


z \ {y — x) °‘dzdwdx, 


(Hi) if y < (—t), then 


wt{y) = —r 


a^i-py-^ 


1—— ft ft fOO 



^1 — Oi 


2paT[l — a) J—tJ\x\Jo 

( w + X 
-i 

2pc 


-Z rn 


w — X 


,2(1-p)' 


-z\ {x — y) °‘dzdwdx. 


Here ra{y) is a density of a positive a-stable random variable with the Laplace transform 


Eg-nZ. ^ g-«“ 

and aV b = max (a, b) 

Proof. From Remark 4.2 in m we get 

P {Y{t) = dy,R{t) = dr) = f f U{dx,dw)vL^,Sa{dy - x,dr + t - w) 

J icGR j wG[0^t] 

where R{t) is a process that counts how long the process Y{t) remains constant from the moment 
t. After applying Eg. (12.31) the above equation reads 

P(F(t) = dy,R{t) =dr) 

= / U{dx,dw) p5dr+t-w{dy - x)vsa.{dr+ t - w) 

JxGR. j iPG[0,i] ^ 

+ (1 - p)6_dr-t+w{dy - x)iysc{dr + t-w) 
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Now we can recover the distribution of Y(t) as a marginal distribution. We calculate the density 
wt{y) applying Tonelli’s theorem to change the order of integration: 


wtiy) = 


a 


r(l-a) 

/ / / 


u{x,w) p6r+t-w{y - x){dr+ t - w) ^ " 
rE[0,oo) ^ 

+ (1 - p)6-dr-t+wiy - x){dr + t - 
t rt 


dwdxdr 


a 


r(l - a) 


[ [ u{x,w) pl(-^,y){x)ll^t-y+x,^){w){y - x) 

J — t ■/ |fc| ^ 


+ (1 - P)'i-iy,oo){x)l(^t+y-x,oo){w)ix “ ?/) 


— l — OL 


dwdx. 


The area of integration in the above integral depends on the relation of y with t - we have three 
cases. The first case is described by the condition \y\ < t. Then 


wtiy) = 


pa 


ry rt 


u{x,w){y — x) ^ °^dwdx 


r(i - a) 
{l-p)a 


\x\\/(t-y+x) 


t rt 


+ 


r(l a) Jy J^x\V{t+y—x) 


u(x, w)(x — y) °^dwdx 


Y y >t then 


wt{y) = 


pa 


ry rt 


u{x,w){y — x) ^ °‘dwdx 


r(l-a) J-t. 


\x\V{t-y+x) 


and finally if y < (—f) then 


wt{y) = 


{l-p)a 


ff 

JV J Ixh 


u{x, w){x — y) °^dwdx 


r(l a') Jy j]^x\\/{t+y—x) 

Substituting Ea. (l2.8p for u{x,w) in all three cases ends the proof. 

In the special case a = ^ we get the following simple expression for wt{y)' 
Corollary 3.1 When a = ^ the PDF wt{y) can be expressed as 

p(i - p)t^ ^ ((^ “ + y)^) +y{it-y)^ + it + y)^) 


wt{y) = 


^ y{t - y) 2 (t + y) 2 ((i - 2p(i -p))t + {i- 2p)y) 


if |y| < t and as 


ifyPt and as 


ify< i-t). 


wt{y) = 


p 

TT 




y (p(y-i)2 + (1 -p)(y+ t)2) 


Wtiy) = 


p-i 


1 

t2 


-y )2 + (1 -p)(-y -t) 2 'j 


TT 


□ 
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Proof. The density Wt{y) is symmetric in a following sense: 


wt{y)p = wt{-y)i-p, 


where wt{y)c denotes the PDF of the process Y{t) with parameter p = c. Thus it is enough to 
calculate wt{y) for y > 0 - we assume now y > 0. There are two cases: y < t and y >t. In the first 
case from Theorem TT\ and Eq. (12.9p we have 


wtiy) 


^ pHi-p) r r* 

2^/^7r JJ^x\\/{t—y+x) 

{p^{w — y) + (1 — p)^{w + y)) (y — x)~^^‘^dwdx 

, p{^-p? r f 

2^/ TT Jy j^x\\/{t+y—x) 

{p^{w — y) + (1 — p)‘^{w + y)) — y)~^/‘^dwdx 


Notice that in the above integrals |x| V {t — y + x) = t — y+x when x > and |x| V {t — y+x) = |x| = 
—X if X < Similarly |x| V (t + y — x) = t + y — x when x < and |x| V (t + y — x) = |x| = x if 
X > Thus Wt{y) can be expressed as a sum of four integrals, each of them having a simple region 
of integration. We can calculate all of them using standard integration techniques. After tedious 
computations we finally get the desired result for wt{y) in case when y G (0, t). In the second case, 
that is y >t, from Theorem 13.11 and Eq. (12.91) we have 


wtiy) = 


P^O^-P? 

2^1'^'k 

r-O rt 


\I I ~ ~ ^^‘^dwdx 

+ j J [p"^{w - y) + {I - p)"^ {w + y)) - x)~^^^dwdxj 


Again, applying standard integration techniques we obtain the desired result for wt{y) when y > ^ 
In a special case a = ^ and p = ^ the PDF wt{y) has the form 

1 (t {{t + y)5 -{t- y)5^ -y{{t- y)^ + {t + y)5)) 


wt{y) = 


27r 


if |y| < t and 


wtiy) = 


TT 


y{t-y) 2 {t + y)h 


1 

t2 

\y\ {\y-t\^ + \y + t\^) 


if \y\ > t. Figure H] presents wi{y) for different values of p. Notice that in all cases we have a 
characteristic sharp peak at y = t and y = —t. However when p —)• 0 the left peak gets bigger and 
the right one diminishes. The opposite situation appears when p —)• 1. 
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Figure 4: For a = 0.5 we compare the densities wi{y) obtained in Corollary 13 .1 1 for different values 
of p (red solid lines) with densities estimated using Monte Carlo methods (blue pluses). 



Figure 5: Plot of PDF wt{y) of the process La{S^^{t)) for different values of t, a = 0.5, p = 0.5. 
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